Excess US Firearm Mortality During the COVID-19 Pandemic Stratified by Intent and Urbanization

This cross-sectional study used time series forecasting to estimate excess firearm mortality in the US during the COVID-19 pandemic.


Hierarchical ensemble time series forecasting
Following the approach of previous literature, [15][16][17] we used time series forecasting to create a counterfactual scenario in which firearm injury followed pre-pandemic trends. We used a standard structure where we defined our training set to be from January 1999 to March 2020 and our testing dataset to include monthly data from April 2020 to December 2021. This structure is described in Figure S1. Using just the training dataset, we developed a time series model reflecting over-time trends in firearm fatalities. Then using this model, we predicted (or forecasted) firearm fatalities monthly for the months of the test dataset (i.e., April 2020-December 2021). For each prediction, we also developed 99% prediction intervals which creates an uncertainty interval around each predicted value. These forecasted values (known as our counterfactual forecasted values) were then compared to the actual measurements observed in April 2020-December 2021.
For modeling we used all test data across the urbanization and intent strata. By placing all strata in the model together, we could also assess trends in the overall firearm fatality rates. We did this by assuming a hierarchical structure with two levels (one level as the overall dataset and another level being each urbanization and intent category).
We applied both autoregressive integrated moving average (ARIMA) and exponential smoothing (ETS) modeling to our training dataset. Briefly, an ARIMA model is specified by the form ARIMA (p,d,q)(P,D,Q)s, where p, d, and q are the autoregressive, differencing, and moving average components of the model, respectively; P, D, and Q are seasonal variants of these components; and s is the order of seasonality. 20 ETS models capture error, trend, and seasonality components of a time series, each of which may be additive or multiplicative in nature; thus, accounting for all combinations, a total of 35 ETS model structures are possible. 21 All model selection was performed using the 'fable' 4 and 'fabletools' 5 RStudio packages, which perform automatic model selection (p, d, q, P, D, Q, and s for ARIMA models; E, T, and S for ETS models) based on minimization of Akaike information criteria (AIC); model selection via these packages have been used in numerous public health forecasting studies, 23-25 including of injury mortality data. 26 To incorporate both ARIMA and ETS components into the models generating our counterfactual forecasts, we used ensemble time series modeling. Ensemble time series modeling incorporates components of multiple time series models into a single model to improve predictive accuracy. The predictive benefits of ensemble (as opposed to single model) forecasting approaches are derived from their inherent ability to capture different patterns present within real-world time series data. 22 In our ensemble model, ARIMA and ETS components were weighted based on the inverse of their variance. This method has the benefit of minimizing the influence of a model component contributing a large error to the forecast and therefore reduces final forecasts' confidence intervals. ARIMA/ETS model parameters were estimated using maximum likelihood method, 4 with fit residuals for the final, aggregate model approximating a normal distribution; Shapiro-Wilks test of fit residuals was non-significant. Ljung-box test was insignificant for the final aggregate model, indicating serial correlation among residuals was random and equal to white noise. 27,28

Model validation
To assess our ensemble model's predictive ability, and therefore its capacity to generate realistic counterfactual trends, we used time series cross validation (TSCV). Cross validation assesses the accuracy of a given model by sub-setting time series data into training and test sets as described above. TSCV does this cross validation procedure on a rolling basis, meaning that we repeat the cross validation process many times using slightly different training and tests sets. 29 Specifically, we first fit ARIMA/ETS ensemble models to five year (January 1999 to December 2004) of pre-pandemic HTS data and forecasted (predicted) 21 months into the future (January 2005 to September 2006); this was the length of pandemic data that was ultimately forecasted. These 21 months of forecasted data (predictions) were then compared statistically to the real data for the same period. Next, the training window was extended by one month (making it February 1999 to January 2005) and 21 months of data (February 2005 to October 2006) were forecasted and again compared to extant data. This process was continued throughout the entirety of the pre-pandemic period. Using this process, we also compared our ensemble ARIMA/ETS model to non-ensemble ARIMA and ETS model.
TSCV-generated forecasts were compared to real data using mean absolute percent error (MAPE). MAPE is an absolute measure of error divided by the corresponding observed value y: MAPE is unit-free, always positive, and may be used to compare forecast accuracy across different datasets. 30 Smaller values of MAPE are preferred.